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ABSTRACT 

The sources of ultra-high energy cosmic rays are not yet known. However, 
the discovery of anisotropic cosmic rays above 57 x 10 18 eV by the Pierre Auger 
Observatory suggests that a direct source detection may soon be possible. The 
near-future prospects for such a measurement are heavily dependent on the flux 
of the brightest source. In this work, we show that the flux of the brightest 
source above 57 x 10 18 eV is expected to comprise 10% or more of the total flux 
if two general conditions are true. The conditions are: 1.) the source objects 
are associated with galaxies other than the Milky Way and its closest neighbors, 
and 2.) the cosmic ray particles are protons or heavy nuclei such as iron and 
the Greisen-Zatsepin-Kuz'min effect is occurring. The Pierre Auger Observatory 
collects approximately 23 events above 57 x 10 18 eV per year. Therefore, it is 
plausible that, over the course of several years, tens of cosmic rays from a single 
source will be detected. 

Subject headings: acceleration of particles 



1. Introduction 

The field of ultra-high energy cosmic rays has long held more questions than answers. 
It has be en over 40 ye ars since the detection of the first cosmic ray with energy greater than 



10 20 eV (ILinsleyi I19631 ) . and we still do not know where these particles come from. However, 



we may be on the verge of rapid progress. 



The detection of anisotropic cosmic rays (CR) above 57 x 10 18 eV (Abraham et al.ll2008 



20071 ) by the newly constructed Pierre Auger Observatory ([Abraham et al.l 12004 ) has lead 
to speculation that CR astronomy may soon be a reality. If this is true, it is expected that 
the brightest source will be the first to be resolved. Since any meaningful astronomical 
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measurement (e.g., sky position, energy spectrum, or morphology) will require the detection 
of many particles from the source object, the near-future prospects for CR astronomy depend 
largely on the flux of the brightest CR source in the sky. Consi der that the Auger Obs ervatory 
detects approximately 23 CR above 57 x 10 18 eV per year ([Abraham et al.l 120071 ). If the 
brightest source contributes 10% of the flux, the Pierre Auger Observatory will detect 23 
events from this source in a decade; the era of CR astronomy will begin. If, on the other 
hand, the flux of the brightest source is far below 1%, the Pierre Auger Observatory may 
not detect more than one CR from this source during its 20-year operational lifetime. 

In this letter, we estimate the flux of the brightest source above 57 x 10 18 eV by postu- 
lating two general conditions. The conditio ns are consistent with a l arge number of leading 
source scenarios; e.g . , rad i o loud galaxies (IFraschetti &: Melial 120081 ). active galactic nuclei 



(IBecker fc Biermannll2009l ; iFarrar &: Gruzinovll2008l ). and gamma ray bursts (jWaxmanlll995 



Wick et al. 12004 ). The first condition is that the source objects are associated with galaxies 



but the sources are relatively rare or rarely active such that many galaxies, such as our 
own and our closest neighbors, do not contain a luminous source. This implies a source 
number density dN/dV < 10~ 2 Mpc -3 . The second condition is that the CR particles 
are protons or heavy nuclei such as iron and the Grei sen-Zatsepin-Kuzmin (GZK) energy 
loss processes ( IGreisenl Il966l ; IZatsepin &: Kuzminl Il966l ) are occurring. Above 10 20 eV, the 
energy loss length for protons and iron is tens of Mpc. This effectively limits the propaga- 
tion distance of the highest energy CR to 100 Mpc, which implies a source number density 
dN/dV > 10" 6 Mpc" 3 . The energy loss lengths of intermediate weight nuclei are less than 
either protons or iron-like heavy nuclei. Therefore, the second condition is conservative, 
i.e. it implies the broadest range of source number densities. This r ange o f so urce densities 



fully e ncompasses the ranges recently estimated by iTakami fc Satol (120091 ) and ICuoco et al. 
5x3). 



2. Simplified Analytical Treatment 

In this work, we only consider the integrated CR flux above a threshold energy E t h = 
57 x 10 18 eV. Because of GZK energy losses, sources beyond a distance D, the so-called 
GZK horizon, will be nearly undetectable in this energy range. The GZK horizon is a strong 
function of the threshold energy, decreasing as E t h increases. By an acc ident of nature, th e 
GZK horizon is approximately the same for both protons and iron nuclei (IHarari et al.ll2006l ). 
Therefore, we will only examine the proton case in detail here, realizing that the iron case 
will have similar results. 



To reproduce the observed energy spectrum we assume that the proton injection spec- 
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tr um (i.e., the n e ar-sou rce spectrum) is a power law dn/dE oc E a with a = 2.6, as suggested 



bvlAUard et al.l (120071) . Our estimation of the GZK horizon then proceeds as follows. As in 
Harari et al. ( 20061 ). we define a GZK attenuation factor A(E t f l ,x) as the fraction of parti- 



cles originally above that still have an energy above E t h after traversing a distance x. 
The attenuation factor is then A(E t h,x) = ((E t h + AE) / 'Eth) 1 ^ ', where AE is the average 
energy loss of a proton traversing a distance x. The energy loss AE is calculated by solving 
dE /dx = —E/\, where A i s the proton energy loss length in the extragalactic medium (e.g., 



sec 



Protheroe fc Johnson! (119961 )). This is the so-called continuous energy loss approxima- 
tion. In Fig. [H we show the results of this calculation for E t h = 57 x 10 18 eV. There is 
a rather sharp breakpoint at approximately 250 Mpc beyond which the attenuation factor 
rapidly decreases exponentially. A reasonable choice for the GZK horizon is then 250 Mpc. 

Let us examine a simplified scenario where there are A" sources of equal luminosity L 
evenly distributed within a sphere of radius D = 250 Mpc centered at earth. We approximate 
the attenuation factor as a step function, with a value of one within the sphere and zero 
outside. The flux received at the detector from a single source inside the sphere is q = uLr~ 2 , 
where u is the detector's sky exposure at the source location, and r is the source distance. 
We assume u> is constant over all sky locations; i.e., the detector sees equally well over the 
entire celestial sphere. 

Because the sources are distributed evenly, the number of sources in a shell of thickness 

dr is 

dN = 3Nr 2 dr/D 3 , (1) 
and the total flux is q to t = $q qdN = 3ujLN / D 2 . We find the median distance to the nearest 

source R by solving J Q dN = 1/2, which gives R = D(2N)~ 1 / 3 . Then the expected flux of 
the nearest source relative to the total flux, is 

Q « uLR- 2 /q tot = (27N/A)- 1 / 3 . (2) 

Since A^ oc pD 3 , where p is the density of sources, Eq. [2] implies Q oc p~ 1 ^ 3 D~ 1 . Because of 
the quadratic form of Eq. []], the median distance is an overestimate of the mean (expected) 
distance. Therefore, Eq. [2] is a lower limit of the expected flux of the brightest source. 
For a source density of 10~ 6 Mpc" 3 , Eq. [2] gives Q = 13%; while for a source density of 
10~ 2 Mpc" 3 , Q = 0.6%. 



3. Numerical Treatment 



In the following numerical treatment, we calculate Q directly via a Monte Carlo al- 
gorithm with several refinements over the analytical treatment: 1.) the GZK attenuation 
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function A is a smooth function of source distance (not a step function), 2.) the sky exposure 
uj is not constant over the entire sky, 3.) the source luminosities are not necessarily all equal, 
and 4.) the distribution of sources is not homogeneous. These refinements tend to increase 
Q above the estimate given in Eq. [2J 

We calculate the flux from each source as q = uALr~ 2 , the attenuation factor A being 
a smooth function of source distance as calculated in the previous section, and the sky 
exposure u j bein g a function of the source declination. We use the uj function developed by 
Sommerd (120011 ) with a detector at latitude —35° (the location of the Southern site of the 
Auger Observatory), and a maximum zenith angle acceptance of 60° (a standard quality cut 
implemented by the Pierre Auger Collaboration). 

We describe the source luminosity distribution with a differential luminosity function 
<p = dN/dL. We consider two scenarios. The first scenario is that is a power law <fi oc L~ 2 
with L ranging over three orders of magnitude. This scenario is motivated by the fact that 
the luminosity function of many astronomical objects (st ars, normal gal a xies, a nd AGN) i s 
a power law over a sim ilar range of luminosities; e.g., see iBasu fc Jonesl (120041 ); IXul (120001 ); 
Mauch fc Sadler! (120071 ). The simple ubiquity of the power law form makes this scenario 
worthy of consideration. The second scenario is the same as in the analytical treatment, 
that all sources have equal luminosity. 

One of our starting assumptions was that the sources are associated with galaxies. The 
distribution of galaxies on Mpc scales is far from homogeneous. Therefore, the number of 
sources in a shell of radius r is not necessarily well approximated by Eq. UJ The largest 
relative deviations from Eq. [1] occur at short distances. For exa mple, the number o f bright 
galaxies (Mb < —19.5) within 100 Mpc is approximately 23,000 (IMarzke et al.lll9981 ). From 
this, Eq. [H predicts that there are 0.023 bright galaxies w ithin 1 Mpc and 2.9 within 5 Mpc. 
However, the observed numbers are 2 and 8, respectively (lTulfyill988l ). To consider this local 
over-density in a simple way, we use a modified radial density function: 



dN/dr 



( 30Nr/D 3 
[ 3Nr 2 /D 3 



: r < 10 Mpc 
: r > 10 Mpc. 



(3) 



We do not try to replicate the directional distribution of galaxies since this has no bearing 
on the results. 

Our numerical treatment calculates Q by averaging the results of many Monte Carlo 
realizations. We form a single Monte Carlo realization as follows. N sources are assigned a 
random sky position, distance, and luminosity. The sky positions are chosen to follow a flat 
distribution (each sky position has equal probability of containing a source), the distances are 
chosen to follow the radial density function given by Eq. [3] out to a maximum distance of 350 
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Mpc, and the luminosities are chosen to follow either a power law distribution (scenario 1) or 
to be all equal (scenario 2). The flux q from each source is computed, then the relative flux 
of the brightest source is computed as Q = max{g 1; . . . , q N }/sum{qi, . . . , q N }. We calculate 
Q by averaging over 1000 realizations. This is repeated for several different source densities. 

We show the results of this analysis in Fig. [2j The error bars represent the 10-90% 
quantiles of each set of realizations. At a source density of 10~ 6 Mpc -3 , the expected value 
of Q for scenario 1 is 40% and the 10-90% quantile range is 19% to 70%. At a source density 
of 10~ 2 Mpc -3 , the expected value of Q is 10% and the 10-90% quantile range is 1.8% to 
26%. 

The numeric results are greater than the analytical treatment at all source densities. 
This is expected since the analytic treatment is a lower limit. From a source density of 
10~ 6 Mpc -3 to 10 -4 Mpc -3 , the numerical results follow the general trend of the analytic 
treatment: decreasing by a factor of approximately 4 as the source density increases by a 
factor of 100. At source densities greater than 10 -4 Mpc -3 , the numerical results flatten 
(i.e., they are nearly independent of source density). This is the source density above which 
we expect the closest sources to be within 10 Mpc and the linear region of Eq. [3] to be 
important. If no local overdensity is assumed, i.e. Eq. [1] is substituted for Eq. [3X then the 
numerical results do not flatten. In this case, the expected value of Q for scenario 1 at a 
source density of 10" 2 Mpc" 3 is 4.0%. 

The relatively weak dependence of Q on source density is caused by two effects which 
somewhat balance each other. The first is that a greater number of sources will tend to 
diminish the relative flux of the brightest source by increasing the background. The second 
is that a greater number of sources will tend to increase the flux of the brightest source since 
there is a greater probability of a source being relatively close or luminous or both. For the 
special case dN/dr oc r, the tendencies exactly balance such that Q is independent of source 
density. 

The weak dependence of Q on source density implies that estimates of source density 
based solely on even t clustering will have a larg e uncertainty. Indee d, the recent CR source 



density estimates by iTakami &: Satol (120091 ) and ICuoco et al.l (120081 ) reflect this 



In general, scenario 1 has a greater Q than scenario 2, however the difference is well 
within the 10-90% quantile range. That is, the results are not sensitive to the details of the 
source luminosity function. Regardless of this, scenario 1 (power-law luminosity function) 
is the more appropriate and the more general consideration. For instance, in scenario 1 the 
brightest sources are spread over a relatively large range of distances, whereas in scenario 
2 the brightest sources are always clustered relatively nearby. If scenario 2 is adopted, a 
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lack of events from a nearby collection of matter can be u s ed to argue against a model with 
relatively high source densities (e.g., see iGorbunov et al.l (120081 )). Such a constraint is the 



result of attributing a difference of magnitude to a difference in type, and is almost certainly 
artificial. 



4. Prospects for CR Astronomy 

We expect that the flux at earth from any particular extragalactic CR source is nearly 
constant over a time span of several years. Even if the particles are accelerated in a rather 
short burst, the CR pulse is broadene d in time becau se of magnetic deflections, with the 



highest energy particles arriving first. IWaxmanl (119951 ) has shown that a time broadening 



of 100 years or more for a 100 Mpc distant source is a reasonable expectation based on our 
present knowledge of extragalactic magnetic fields. If source persistence is indeed due to 
extragalactic magnetic fields, the apparent source density and apparent source luminosity 
are functions of these intervening magnetic fields, i.e. they are functions of distance. In 
this case, the local apparent source density should be used with Fig. [2] since we are only 
concerned with the flux of the brightest source which is typically one of the nearest. 

If we assume the flux from the brightest source is steady and Q is 10%, then we expect 
the Pierre Auger Observatory to acquire 2.3 events from this source per year. The angular 
size of the source is dependent on the magnetic rigidity of the particle and the intervening 
magnetic fields, both of which are not well constrained. However, even if the source covers a 
20° x 20° area (approximately 1% of the sky observed by the Pierre Auger Observatory), the 
number of events from this region will be ten times the expected amount from an isotropic 
background. In this respect, a source detection is expected to occur in a matter of years. 

In this work, we have concentrated on the energy range above 57 x 10 18 eV. However, if 
the extragalactic source paradigm is correct, CR astronomy is possible down to the energy 
where the diffuse galactic CR become domi nant over the extragalactic CR (most likely 



somewhere below 10 19 eV (lAUard et al.l 120071 )). At lower energies, Q will be less (because 



as the energy threshold E t h decreases the GZK horizon D increases, and Q is inversely 
proportional to D) and the angular size of the sources will be larger because of decreasing 
magnetic rigidity. 



The author thanks his fellow Pierre Auger collaborators, and acknowledges the support 
of the National Science Foundation (award number 0838088) and the Michigan Space Grant 
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Fig. 1. — The GZK attenuation factor of a proton source observed in an energy range E > 
57x 10 18 eV. We assume that the near-source energy spectrum is a power law dn/dE oc E~ 2S . 




Fig. 2. — The expected flux, relative to the total flux, of the brightest cosmic ray source 
above 57 x 10 18 eV as a function of source number density. The points are slightly offset on 
the x-axis for clarity. The analytic approximation is Eq. [2j 



